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We develop an elementary mean field approach for fully asymmetric kinetic Ising models, which 
can be applied to a single instance of the problem. In the case of the asymmetric SK model this 
method gives the exact values of the local magnetizations and the exact relation between equal- 
time and time-delayed correlations. It can also be used to solve efficiently the inverse problem, 
i.e. determine the couplings and local fields from a set of patterns, also in cases where the fields 
and couplings are time-dependent. This approach generalizes some recent attempts to solve this 
dynamical inference problem, which were valid in the limit of weak coupling. It provides the exact 
solution to the problem also in strongly coupled problems. This mean field inference can also be 
used as an efficient approximate method to infer the couplings and fields in problems which are not 
infinite range, for instance in diluted asymmetric spin glasses. 
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Inference problems are as old as scientific modelling: given data, how can we reconstruct a model which accounts for 
it, and find the parameters of the model? This is particularly difficult when data is obtained from networks of many 
interacting components. The fast development of high-throughput technologies in various fields of biology, ranging 
from gene regulation to protein interaction and neural activity, is generating a lot of data, which is challenging our 
ability to infer the structure and parameters of the underlying networks. 

This 'network reconstruction' problem is typically an inverse problem which has motivated a lot of activity in 
machine learning and in statistical physics [^]-^ p[-p^ p^|-pH p3| , p6| , p8| , pof . Until recently the main efforts have 
been dedicated to reconstructing equilibrium Boltzmann-Gibbs distributions. In the so-called inverse Ising model, 
one typically assumes to have data in the form of some configurations, which we shall call 'patterns', of a TV-spin Ising 
system drawn from the Boltzmann-Gibbs distribution with an energy function including one-body (local magnetic 
fields) and two-body (exchange couplings) terms. The problem is to reconstruct the local fields and the exchange 
couplings (collectively denoted below as 'couplings') from the data. This problem has been actively studied in 
recent years, in particular in the context of neural network inference based on multielectrode recordings in retinas 
[Q, |2^. The standard solution of this problem, known as the Boltzmann machine, computes, for some given 
couplings, the local magnetizations and the two-spin correlations, and compares them to the empirical estimates of 
magnetizations and correlations found from the patterns [Q, ^]. The couplings are then iteratively adjusted in order 
to decrease the distance between the empirical magnetizations/correlations and the ones computed from the model. 
A Bayesian formulation shows that the problem of finding the couplings is actually convex, so that this iterative 
procedure is guaranteed to converge to the correct couplings, provided that the number of patterns is large enough 
to allow for a good estimate of correlations. The drawback of this method is that the reliable computation of the 
magnetizations/correlations, given the couplings, which is done using a Monte Carlo procedure, is extremely time 
consuming. Therefore this exact approach is limited to systems with a small number of spins. Most of the recent 
works on this issue have developed approximate methods to infer the couplings. Among the most studied approaches 
are the naive mean field method 13, 2^, the TAP approach a method based on a small magnetization 



expansion and a message-passing method called susceptibility propagation ||TT|, ^ Another approach which 
has been developed is that of linear relaxation of the inference problem [g^. The inverse-Potts problem is a version 
of this same problem, with variables taking q possible states. The case g = 20 is relevant for inferring interaction 
in protein pairs from data on co-evolution of these pairs, and its solution by susceptibility propagation has given 
an accurate prediction of inter-protein residue contacts p8|. Another case which has received some attention is the 
problem of reconstruction in Boolean networks (see e.g.g^and references therein). 

However, in many applications to biological systems, in particular the ones concerning neural activity and gene 
expression network, the assumption that the patterns are generated by an underlying equilibrium Boltzmann-Gibbs 
measure is not well founded. Couplings are typically asymmetric, and they may vary in time, so there is no equilibrium 
measure. This has prompted the recent study of inference in purely kinetic models without an equilibrium measure 
iQ, 1^, 2^. A benchmark on this dynamic inference problem is the inverse asymmetric kinetic Ising model. The 
framework is the same as the equilibrium one: one tries to infer the parameters of the dynamical evolution equation 
of an Ising spin systems, given a set of patterns generated by this evolution. The recent works ^, ^0|, |^, ^ have 
studied the performance of two mean-field methods on this problem, the naive mean field (nMF) and a weak-coupling 
expansion which they denote as TAP method. They have shown that, in the case of the fully asymmetric infinite 
range spin glass problem, the inference problem can be solved by these methods in the case where the spins are weakly 
coupled. In the present work we present a (non-naive!) mean field approach which solves the problem at all values of 
the couplings (and reduces to their TAP approach at weak coupling). 

The kinetic Ising model which we shall study is the same as the one of N Ising spins Si evolve in discrete time, 
with a synchronous parallel dynamics. Given the configuration of spins at time t—1, s{t—l) = {si(i— 1), . . . , spfit—l)}, 
the spins Si{t) are independent random variables drawn from the distribution: 

^ 1 

^'""i'"-'»-n ,,„i.(^t.w) '"""'""' <" 

where 

h^{t) ^ H^{t-l)+Y^ .hj{t - l)Sj{t - 1) (2) 
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Note that both the local external fields Hi{t) and the exchange couplings Jij{t) may depend on time. Here we are 
interested in a fully asymmetric model. We generate the Jy by an asymmetric version I, § of the infinite - 
range Sherrington-Kirkpatrick spin glass model |2^ , in which for each directed edge (ij) the coupling is a gaussian 
random variable with variance 1/7V. Notice that Jij and Jji are independent random variables. We do not include 
self-interactions (we take Ju = 0), although this could be done without changing the results. As initial conditions we 
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take Si{t = 0) = ±1 with probability 1/2. Our method also applies to the case of asynchronous dynamics, studied 
in with the TAP approach, but to keep the presentation simple we shall study only the case of the synchronous 
parallel dynamics in this letter. 

We first derive the mean- field equations for the magnetizations mi{t) = {si{t)). Because the couplings are asym- 
metric, JijJji = 0(1/ y/N), therefore the Onsager reaction term is not present in this problem. This makes the 
derivation of our equations, which we shall denote just 'mean-field' equations, particularly easy. The approximate 
equations used in |29| , originally derived in fl^ , have been obtained by a second order expansion in the couplings. 
When this expansion is applied to the symmetric problem it gives back the TAP equations p7[ with their Onsager 
reaction term. In the present case of asymmetric coupling, it keeps the correction of order Jij Jij . We shall keep 
for these second-order-expanded equations the name 'TAP'-equations, as used by jl^, |2^, [29| ]. 

The local field on spin i due to the other spins, Jij(t — l)sj(t — 1), is the sum of a large number of terms. 
Therefore it has a gaussian distribution with mean 

j 

and variance 

A,it-l)=J24{l-m,it-m (4) 

j 

(in order to derive this last formula, one must use the fact that the typical connected correlation (sjSk) — mjirik is 
typically of order I/^/N; this will be checked self-consistently below). Using this property and the definition of the 
dynamics (|l|), one obtains the magnetization of spin i at time t: 

m,{t) ^ j Dx tanh [/3 (H,{t - I) + g,{t - 1) + x^ A,{t - 1)'^'^ , (5) 

where Dx — -j^^ ^ is the probability density for a Gaussian variable x with zero mean and variance unity. 

Equations (^|^,^ are our mean field (MF) equations for this problem, valid on a given instance. Similar dynamical 
equations have been obtained in the study of the sample-averaged order parameter in asymmetric neural networks]^, |^ 
and spin glasses]^. They can be iterated starting from some initial condition (in our case TOi(O) = 0) in order to get 
all the magnetizations mi{t) at any time. They rely only on the central limit theorem and they are exact in the large 
iV limit, for any set of couplings and external fields, even if they are time-dependent. These differ from the 'TAP' 
equations of ||, |l[ |§ which can be written in our notation: 

m,{t) ^ t&nh [13 H,{t-l) + /3g,{t~l)-m,{t) (3^ A,{t~l)] , (6) 

and from the naive mean field (nMF) equations: 

m,(t) = tanh [/3 {H,{t - 1) + g,{t - 1))] . (7) 

The nMF equations and the 'TAP' equations actually give the same result as our exact MF equations, when expanded 
in powers of A^, respectively to order A? and Aj, but they differ at order Af . The fact that 'TAP' equations agree 
with the exact MF to second order in a weak coupling expansion is consistent with their derivation through second 
order Plefka-type expansion The correctness of the MF equations (|^,^,^ can be easily checked numerically as 
shown in the left panels of Fig^. 

We now turn to the computation of correlations. We shall establish the mean field relation between the time-delayed 
and the equal-time correlation matrices: 

D,,{t)^{5s,{t + l)5s,{t)) (8) 
C,,{t) = {Ss,{t)Ss,{t)) , (9) 

where we define 6si{t) as the fluctuation of the magnetization: Ssi(t) = Si{t) — {si(t)) . 

We start by writing Jij{t)sj{t) — gi{t) + 5gi{t), where Sgi{t) is gaussian distributed with mean and variance 
Ai{t). Now, by definition of Dij we have 

D,j{t) = {s,it) tanh [/? iH,{t) + g,{i) + 5g,{t))]) - (s,(i)>(tanh [/3 {H,{t) + g,{t) + 5g^m) (10) 
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Hereafter in order to keep notations simple in the derivation of the relation between D{t) and C{t) we work at a fixed 
time t and we thus drop the explicit time indices: all time indices in this derivation are equal to t (e.g. = Jijit), 
Ssi = Ssi{t), Qi = gi{t) etc.) We get: 



JjkDik = {{gj + 6gj) tanh [/S {Hi + gi + 6gi)]) - gj{tanh [P {Hi + gi + %)] 

k 

= {6gj tanh [(3 {H, + g, + 5g^)]) 



(11) 



In order to evaluate the average we need the joint distribution of Sgi and Sgj. The crucial point to keep in mind 
is that, as the couplings are of order each matrix element of C and D is also small, of order 1/^/N. Their 

covariance is therefore small: 



{'>9iS9j) = (Y^ Jtk {sk - (sfc)) ^ Jji {si - (s;))) 

k I 
k,l 



(12) 
(13) 



where e is typically of order l/yiV. So the joint distribution of a; = Sgi and y = Sgj takes the form, in the large N 
limit (omitting terms of order e^): 



P{x,y) 



1 



: CXp 



y 

27ryA^ """" V 2A, 2Aj 
Using the small s expansion of eq. (^4|) we can rewrite eq. (|ll|) as 

k 



_xy_ 



2 A,- 2A, 



dx 



xy^ tanh [/3 (iJi + gt + x)] 



J ^/2¥A: 

Combining eq. ( |l^ ) and eq. (^5|) we get: 

(i^J^),„. = (JCJ^)^,./3 



: exp (1 - tanh^ [/? {H, + g, + x)]) 



dx 



:e ^ (l - tanh^ /3 (iJj + + x)) , 



(14) 

(15) 
(16) 

(17) 



which gives the explicit mean-field relation between C and D. Putting back the time indices, we obtain the final 
result in matrix form: 



D{t) = A{t) J{t) C{t) 
where A{t) is a diagonal matrix: Aij{t) — ai{t)6ij, with: 



(18) 



(19) 



The final result (|l8| ) takes exactly the same form as the one found with the naive mean field equation and with 
the 'TAP' approach. The predictions of all three methods, uMF, 'TAP' and our MF method is always D{t) = 
A{t) J{t) C{t), with a diagonal matrix A{t) which differs in each case. As shown in ||2^, the nMF approximation 
gives: 



a,{t) ^13 f Dx [l - tanh^ /3 (^H,{t) + g,{t) + x^fK^f)) 



a^^^{t)^P[\-m,{t + lf 



the 'TAP' approximation gives: 



al''^{t + l)^ P\l~m,{t+IY 



and our mean field prediction is the one given in (pj|). 



(l-m,(t + l)2)/?2^J^(l-mfe(t)2) 



(20) 



(21) 
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We claim that, as in the case of the magnetizations, our mean field equations connecting D to C are exact in the 
asymmetric SK model, in the large N limit. This statement can be checked numerically by comparing (AJC)ij with 
the experimental values of Dij found by monte carlo simulations, as shown in Fig.|l|. 

These results on the mean field relation between C and D can be used for the inverse problem. Given P 'patterns', 
which are time sequences of length t generated from the distribution (^, one can estimate for each t — 1, . . . , i, the 
magnetizations mi{T), the equal time correlations Cijlr) and the time-delayed correlations Dij{T). The problem is 
to infer from these data the values of the couplings Jij (r) and of the local fields Hi (r) . Without loss of generality, 
we can use /3 = 1 as it is absorbed in the strength of couplings and fields that we want to infer. We shall solve this 
problem using the mean field equations. 

The problems corresponding to different times and sites decouple. So let us consider a fixed value of i and r, and 
infer the Jijir) for j — I, . . . , N , and Hiir). To lighten notation we drop the explicit indices r and i, and we denote 
H = Hi{T), rrij = mj{T), m = TOi(T + 1), g ~ giir), A = Ai(T), a = ai{T). Following one can obtain J by 
inverting the relation (|l^). The first step is to invert the empirical C matrix and compute: 

6,=5]Afe(T)C^/(r) . (22) 



If one knows the number a = ai{T) one can then infer the couplings from (|18|): 

MT)^h,/a. (23) 
Let us now see how a can be computed. The mean field equation (||) for the magnetization reads: 

J Dx tanh H + g + xVA . (24) 

The equation ( |l^ for a is 

J Dx(l- tanh^ H + g + xVA ) . (25) 
The link between a and A is obtained from (P, which reads: 

i 

To solve this system of equations, we propose the following iterative procedure. Using the empirical correlations 
and magnetizations estimated from the patterns, we first compute from ( p2| ) the {bj}, j € {l,...,iV}, and 7 = 

E,&?(i--I)- 

Then we use the following mapping to find A. 



Start from a given value of A. 

Using the empirical value of m and the value of A, compute H + g hy inverting (p^). The right-hand side of 
this equation is an increasing function oi H + g so this inversion is easy. 

• Using H + g and A, compute a using ( p5| ) 

• Compute the new value of A, called A, using (p6|). 

It is worth pointing out that in the thermodynamic limit, iV — > cxd, the value of A becomes independent of i. So, if 
the system under consideration is large enough, the above iteration could be perfomed only once in order to reduce 
computation time. 

This procedure defines a mapping from A to A = /(A), and we want to find a fixed point of this mapping. It turns 
out that a simple iterative procedure, starting from an arbitrary Aq (for instance Aq = 1) and using A„_|_i = /(A„), 
usualUy converges. More precisely, it can be shown that /(O) — )^^_^2y and that the asymptotic form for the slope 

of / for A ^ 1 is /' ^ ^7exp(M^)A„, where m is such that m = erf('u/-\/2). We have found numerically that when the 
number of patterns is large enough the slope verifies: df /dA g]0, 1[. Therefore the mapping converges exponentially 
fast to the unique fixed point. This method therefore works when the number of patterns per spin P/N is large 
enough. In the double limit P, — 00 and P/N large enough the above procedure thus allows to get the exact result 
for A; and therefore to find the couplings Jijir) = bj/a. Once the couplings have been found, one can easily compute 
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g = J2j Jij{'T)^j{''')j s-nd therefore get the local field H{t). The number of operations needed for the full inference of 
the couplings and fields is dominated by the inversion of the correlation matrix C, a time which is typically at most of 
order N^. If the number of patterns is too small, it may happen that there is no solution to the fixed point equation 
/(A) = A. Then one can decide to use A ~ /(O), which is nothing but the nMF estimate for ai(r). 

We have tested our mean field inference method on the asymmetric SK problem, where the couplings Jy are time- 
independent, gaussian distributed with variance and the fields are time independent, uniformly distributed on 
[—/?,/?] . Fig.(||) shows a scatter plot of the result on one given instance at /3 = .4 and /3 = 1.4, and compares it to 
the inference method of using nMF and 'TAP' (the 'TAP' inference is limited to small values of (3: at large (3 it 
fails). Figs. || and ^ show a statistical analysis of the performance of MF inference. It accurately infers the couplings 
and fields even in the strong coupling regime. 

The method that we propose is exact and allows for a very precise inference of the couplings when applied to the 
fully asymmetric SK spin glass, at any temperature, if the number of patterns is large enough. At the same time, it 
is an easy and versatile method which can be used as an approximate inference method when the number of patterns 
is not very large (although one should at least have P > in order for C to be invertible), or when the underlying 
model is not of the SK type. As an example showing the possible use of the method, we have applied it to a sample 
where the matrix is sparse, generated as follows. We first generate a regular graph with 200 vertices and degree 
6 on each vertex. For each edge ij of this graph we choose randomly with probability 1/2 an orientation, say i ^ j. 
Then we take Jji = and Jij is drawn randomly from the probability density (|x|/2)e~^ All the other couplings 
corresponding to pairs of sites kl which are not in the graph are set to 0. One then iterates the dynamics (|l|) 100000 
times at /3 = .6, and uses this data to reconstruct the couplings. Fig. || shows the resulting couplings as a scatter 
plot. The topology of the underlying interaction graph can be reconstructed basically exactly, both by nMF and MF 
by using a threshold, deciding that all reconstructed couplings with \Jij\ < .04 vanish. The non-zero couplings are 
found accurately by the MF inference method. 

To summarize, we have introduced a simple mean field method which can be applied on a single instance of a 
dynamical fully asymmetric Ising model. In the case of the asymmetric SK model this MF method gives the exact 
values of the local magnetizations and the exact relation between equal-time and time-delayed correlations. This 
method can be used to solve efficiently the inverse problem, i.e. determine the couplings and local fields from a set 
of patterns. Again this inference method is exact in the limit of large sizes and large number of patterns, in the 
asymmetric SK case. It can also be used in cases where the underlying model is different, for instance for diluted 
models. This could be quite useful for many applications. 

We thank Lenka Zdeborova for useful discussions. This work has been supported in part by the EC grant 'STAMINA', 
No 265496. 
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FIG. 1. Magnetizations (left column) and correlations (right column) obtained by MF (blue), 'TAP' (green) and nMF (red). 
One A'" = 200 spin model is simulated 500000 times for 31 time steps. The three plots in each column correspond to inverse 
temperature /3 = .3, .8 and 1.4 (from top to bottom). In the left column, the magnetizations predicted by each method for all 
time steps are plotted versus the experimental ones found by monte carlo simulation. For the plots of the right column, the 
correlation matrices C and D are obtained at t = 30. The scatter plot shows for each pair ij, the value of Dij in ordinate, and 
the value of {AJC)ij in abscissa. The three methods differ in their predictions for A. At high temperature, /3 — .3, all methods 
are good for both the magnetizations and correlations; the MF and 'TAP' methods nearly coincide and are slightly better than 
nMF. At larger and larger /3, the 'TAP' correction to naive mean field overshoots, and only the MF results is correct. The data 
supports the statement that MF is exact at all temperatures, while nMF and 'TAP' are only high temperature approximations. 



FIG. 2. Left: The infered couplings found by MF (blue), 'TAP' (green) and nMF (red) plotted versus the real ones for a 
= 100 model, given P = 1000000 patterns generated at inverse temperature /3 = 0.4. Right: The same for (3 = 1.4 (MF 
(blue) and nMF (red), 'TAP' is not shown as it fails at this high /3) 





FIG. 3. Mean square error of the infered couplings (J|"'"'^'^ — J^j^^)^ obtained by MF inference (blue), 'TAP' (green) and nMF 
(red) versus the number of patterns used to estimate the correlations, for a system of size A*' — 40, where the patterns were 
generated from a gaussian distribution with a root-mean-square /3 = 0.2 (left) and — 0.6 (right). The curves are averages 
performed over 20 realizations of the couplings and fields. Notice that the 'TAP' method is absent in the right figure because 
it fails to provide results at strong coupling. 
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FIG. 4. Mean square error of the reconstructed couplings versus /3, averaged over 10 systems with 100 spins, using the three 
inference methods nMF, 'TAP' and MF, with a number of patterns P = 10000, P = 100000 and P = 1000000. All three 
methods agree at small (3. The nMF error can increase by several orders of magnitude at large /3. The 'TAP' method fails to 
provide results above /3 ~ 0.4. The MF inference method gives good results in the whole range of f5. 
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FIG. 5. Mean field inference of a finite connectivity model. The couplings found by MF (blue) and nMF (red) are plotted versus 
the real couplings used to generate the data for a N — 200 model on an asymmetric random regular graph with connectivity 
c — 6 (average in-degree=3) , given P = 100000 patterns generated at inverse temperature /3 — .6. 



